Patient Enrollment Filtering: Claims data were filtered to include only patients with enrollment periods exceeding 12 months. This criterion was set to ensure that the data represented a sufficiently long period to analyze treatment patterns and outcomes accurately.
Linking rx_information: The rx_information data was linked to the claims data using National Drug Codes (NDCs). This linkage was essential for integrating prescription information with the claims data.
Identifying Buprenorphine OUD Cohort: The buprenorphine opioid use disorder (OUD) cohort was identified using the Universal Service Code (usc_cd) of 78430, which corresponds to ‘Drug Dependence’.
Chunk Processing Function: A custom function was written to process each chunk of data separately. This approach ensured that the dataset was handled efficiently and effectively, especially given the potentially large size of the claims data. The processed chunks were then combined into a single dataset.
Duplicate Removal: Duplicate entries were meticulously removed to prevent redundancy and ensure the accuracy of the dataset.
Joining with Demographics: The resulting data frame was joined with demographic information to create the final dataset, named ‘claims_demo’. This final dataset integrated all necessary data elements, including claims, prescriptions, and patient demographics, providing a comprehensive foundation for subsequent analyses.
The final cohort now has 3358850 claims total and 187790 individual patients.
Linking Demographic Data with Claims: The demographic data was linked to the claims data via patient IDs. This step ensured that each claim could be associated with the corresponding patient’s demographic information.
Linking RX Information with Claims: The rx_information data was linked to the claims data using the NDCs. This linkage was essential for associating specific medications with patient claims.
‘claims_demo’ is the final dataset.
library(dplyr)
library(data.table)
library(ggplot2)
library(plotly)
library(usmap)
library(gridExtra)
Attaching package: ‘gridExtra’
The following object is masked from ‘package:randomForest’:
combine
The following object is masked from ‘package:dplyr’:
combine
filter_cohort<-function(df) {
class(df$ndc) <- "character"
class(rx_information$ndc) <- "character"
claims_rx <- df %>%
left_join(.,rx_information, by = "ndc")
keywords <- c("BUNAVAIL", "BUPRENORPHINE/NALOX", "BUPRENORPHINE HCL", "BUPRENORPHINE HCL/NALOXON", "BUPRENORPHINE HYDROCHLORI", "PROBUPHINE IMPLANT KIT", "SUBLOCADE", "SUBOXONE", "SUBUTEX", "ZUBSOLV")
pattern<-paste(keywords,collapse="|")
claims_rxb <-claims_rx %>%
filter(grepl(pattern, product_name, ignore.case=TRUE) | grepl(pattern, generic_name, ignore.case=TRUE))
claims_rxb_12months <- claims_rxb %>%
filter(.$pat_id %in% patients_full_enroll)
duplicates<-apply(claims_rxb_12months,1,function(row) length(unique(row)==1))
claims_rxb_12months<-claims_rxb_12months[!duplicates, ]
claims_rxb_12mo<-claims_rxb_12months[claims_rxb_12months$usc_cd == "78340", ]
fwrite(claims_rxb_12mo, "df2_filtered")
claims_rxb_12mo2<-claims_rxb_12mo %>% select(pat_id, dayssup, quan, from_dt, to_dt, ndc, proc_cde, diagprc_ind, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8) %>% arrange(pat_id)
fwrite(claims_rxb_12mo2, "df2_filtered_selected")
return(head(claims_rxb_12mo2))
}
filter_cohort(d2)
missingness <- colSums(is.na(claims_demo[,c(1:20)])) / nrow(claims_demo[,c(1:20)]) * 100
missingness_df <- data.frame(
Column = names(missingness),
MissingPercentage = missingness
) %>%
arrange(desc(MissingPercentage))
p <- plot_ly(
data = missingness_df,
x = ~MissingPercentage,
y = ~reorder(Column, MissingPercentage),
type = 'bar',
orientation = 'h',
text = ~paste("Missing: ", round(MissingPercentage, 2), "%"),
hoverinfo = 'text',
marker = list(color = 'steelblue')
) %>%
layout(
title = 'Percentage of Missing Values by Column',
xaxis = list(title = 'Missing Percentage (%)'),
yaxis = list(title = 'Column', automargin = TRUE, tickfont = list(size = 10)),
margin = list(l = 200, r = 50, t = 50, b = 50) # Adjust margins to fit labels
)
p
missingness <- colSums(is.na(claims_demo[,-c(1:20)])) / nrow(claims_demo[,-c(1:20)]) * 100
missingness_df <- data.frame(
Column = names(missingness),
MissingPercentage = missingness
) %>%
arrange(desc(MissingPercentage))
p <- plot_ly(
data = missingness_df,
x = ~MissingPercentage,
y = ~reorder(Column, MissingPercentage),
type = 'bar',
orientation = 'h',
text = ~paste("Missing: ", round(MissingPercentage, 2), "%"),
hoverinfo = 'text',
marker = list(color = 'steelblue')
) %>%
layout(
title = 'Percentage of Missing Values by Column',
xaxis = list(title = 'Missing Percentage (%)'),
yaxis = list(title = 'Column', automargin = TRUE, tickfont = list(size = 10)),
margin = list(l = 200, r = 50, t = 50, b = 50) # Adjust margins to fit labels
)
p
Since diag9-diag12 are missing 100%, we could remove those columns
claims_demo<- claims_demo %>%
select(-diag9, -diag10, -diag11, -diag12)
Since diag9-diag12 are missing 100%, we could remove those columns
Group by patient ID and dates claimed
claims_demo <- claims_demo %>% arrange(pat_id, from_dt) %>% group_by(pat_id)
claims_demo <- claims_demo %>%
mutate(year = year(ymd(from_dt)))
head(claims_demo,100)
Summaries of Days prescribed and quantity prescribed
table(claims_demo$dayssup)
-90 -60 -35 -33 -32 -30 -29 -28
2 4 1 1 1 247 2 24
-27 -26 -25 -24 -23 -22 -21 -20
1 3 1 1 1 2 6 10
-18 -16 -15 -14 -12 -11 -10 -9
2 5 26 30 1 1 21 6
-8 -7 -6 -5 -4 -3 -2 -1
16 26 19 129 18 17 33 28
0 1 2 3 4 5 6 7
7569 35993 43963 45592 40734 59788 40755 392620
8 9 10 11 12 13 14 15
63300 15387 61231 12231 19830 15489 303467 140283
16 17 18 19 20 21 22 23
17614 7371 9899 4454 35576 56336 18637 14852
24 25 26 27 28 29 30 31
12068 26315 14936 14248 391483 16661 1391960 2739
32 33 34 35 36 37 38 39
2205 1121 1818 1670 350 238 141 34
40 41 42 43 44 45 46 47
714 17 331 29 34 721 26 12
48 49 50 51 52 53 54 55
45 18 170 11 14 17 14 14
56 57 58 59 60 61 62 63
245 6 9 4 1703 9 3 12
64 65 66 67 68 69 70 71
20 8 3 4 6 1 30 2
72 73 75 76 77 78 80 82
7 3 54 13 3 3 49 1
83 84 85 86 87 88 89 90
5 73 13 1 1 13 3 3784
92 93 95 96 99 100 105 112
1 5 1 2 1 11 1 3
118 120 128 140 141 150 180 240
1 7 1 1 1 5 18 2
257 300 301 320 328 500 999
1 7 2 1 1 2 3
dayssup_freq <- table(claims_demo$dayssup)
df_dayssup <- as.data.frame(dayssup_freq)
colnames(df_dayssup) <- c("dayssup", "frequency")
n <- nrow(df_dayssup)
part_size <- ceiling(n / 4)
df_dayssup_part1 <- df_dayssup[1:part_size, ]
df_dayssup_part2 <- df_dayssup[(part_size + 1):(2 * part_size), ]
df_dayssup_part3 <- df_dayssup[(2 * part_size + 1):(3 * part_size), ]
df_dayssup_part4 <- df_dayssup[(3 * part_size + 1):n, ]
plot1 <- ggplot(df_dayssup_part1, aes(x = as.factor(dayssup), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 1)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot2 <- ggplot(df_dayssup_part2, aes(x = as.factor(dayssup), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 2)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot3 <- ggplot(df_dayssup_part3, aes(x = as.factor(dayssup), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 3)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot4 <- ggplot(df_dayssup_part4, aes(x = as.factor(dayssup), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 4)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(plot1)
print(plot2)
print(plot3)
print(plot4)
table(claims_demo$quan)
-90 -60 -58 -56 -28 -26 -12 -7 0
2 1 1 1 1 1 1 1 54831
0.5 0.7 0.9 1 1.1 1.111 1.3 1.5 2
2442 1 2 11794 3 3 2 2707 25628
2.5 3 3.5 4 4.5 5 5.5 5.6 6
1 22863 3 31202 1 26540 3 1 32881
6.5 6.99 7 7.03 7.5 8 9 9.5 9.9
1 1 86479 1 4 31446 15973 1 1
10 10.5 11 12 12.5 13 13.5 14 15
56783 9 45648 30625 4 11793 1 254479 87155
16 16.04 17 17.5 18 19 20 21 22
23352 1 7501 2 34466 4403 38616 106244 14743
22.5 23 23.8 24 24.5 25 26 27 27.5
5 16811 1 16383 1 17791 8446 6888 1
28 29 29.7 30 31 32 33 34 35
189114 3517 1 345997 2176 12683 3809 4324 29121
36 37 37.5 38 39 40 41 42 43
7289 4231 1 12975 3054 15890 1593 102241 1562
44 45 46 47 48 49 50 51 52
4142 158437 3579 979 4478 9130 13049 1469 4969
52.5 53 54 55 56 57 58 59 60
4 7944 5445 3614 155650 1660 6326 649 509950
61 62 63 64 65 66 67 68 69
312 1463 9027 1249 2593 1517 753 3054 884
70 71 72 73 74 75 76 77 78
26972 241 1439 760 659 63298 543 1556 1256
79 80 81 82 83 84 85 86 87
209 3551 1063 609 1386 47881 1455 621 1660
88 89 90 91 92 93 94 95 96
737 114 271522 217 165 477 86 278 292
97 98 99 100 101 102 103 104 105
97 1543 78 3707 22 1566 16 131 6096
106 107 108 109 110 111 112 113 114
77 31 239 19 406 42 5471 76 92
115 116 117 118 119 120 121 122 123
182 260 23 48 31 56357 10 28 31
124 125 126 127 128 129 130 131 132
120 65 254 17 60 11 231 2 38
133 134 135 136 137 138 139 140 142
2 34 967 56 3 23 2 489 10
143 144 145 146 147 148 149 150 151
4 71 42 8 29 12 1 6053 1
152 153 154 155 156 157 158 159 160
7 5 37 10 11 3 4 5 64
162 163 164 165 166 167 168 170 172
12 1 9 113 7 1 399 31 10
173 174 175 176 177 178 179 180 181
1 14 7 5 4 3 1 7300 1
182 184 185 186 187 189 190 192 195
11 3 5 7 1 5 15 11 36
196 198 200 201 203 204 205 206 208
267 1 108 1 7 11 2 1 3
209 210 211 212 215 216 217 218 220
6 726 1 1 2 5 1 1 38
221 224 225 228 230 231 232 240 243
1 134 129 1 9 1 13 1757 2
244 245 246 250 252 255 256 257 260
2 1 1 8 17 3 1 1 8
261 264 265 268 269 270 275 280 290
1 1 1 1 1 1267 4 37 11
294 300 304 306 315 320 322 324 325
2 321 1 1 7 3 1 2 2
330 333 336 340 341 348 350 354 360
25 1 15 1 1 3 4 4 577
370 375 380 384 390 400 405 410 420
3 3 22 1 20 4 6 1 29
440 442 448 450 464 476 480 500 504
1 1 1 82 2 1 60 1 1
510 520 540 543 550 560 570 580 600
3 3 44 2 1 23 1 3 162
620 630 660 680 690 700 720 750 780
1 2 1 1 1 3 6 13 3
800 810 840 870 890 900 930 960 980
2 14 12 1 1 62 1 1 5
1000 1050 1060 1080 1100 1120 1200 1230 1800
3 5 2 1 1 5 20 1 2
2000 3000 5000 5829 7000 8000 9000 11000 12000
10 10 8 1 12 4 2 6 10
13000 14000 16000 21000 22000 23000 24000 25000 26000
2 5 6 20 9 12 1 18 1
27000 32000 33000 35000 38000 40000 41000 42000 43000
2 10 4 2 4 6 2 13 1
44000 45000 46000 50000 53000 56000 62000 63000 66000
2 63 3 2 6 1 1 3 2
270000
3
quantities_freq <- table(claims_demo$quan)
df_quantities <- as.data.frame(quantities_freq)
colnames(df_quantities) <- c("quantities", "frequency")
n <- nrow(df_quantities)
part_size <- ceiling(n / 8)
df_quantities_part1 <- df_quantities[1:part_size, ]
df_quantities_part2 <- df_quantities[(part_size + 1):(2 * part_size), ]
df_quantities_part3 <- df_quantities[(2 * part_size + 1):(3 * part_size), ]
df_quantities_part4 <- df_quantities[(3 * part_size + 1):(4 * part_size), ]
df_quantities_part5 <- df_quantities[(4 * part_size + 1):(5 * part_size), ]
df_quantities_part6 <- df_quantities[(5 * part_size + 1):(6 * part_size), ]
df_quantities_part7 <- df_quantities[(6 * part_size + 1):(7 * part_size), ]
df_quantities_part8 <- df_quantities[(7 * part_size + 1):n, ]
plot1 <- ggplot(df_quantities_part1, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 1)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot2 <- ggplot(df_quantities_part2, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 2)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot3 <- ggplot(df_quantities_part3, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 3)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot4 <- ggplot(df_quantities_part4, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 4)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot5 <- ggplot(df_quantities_part5, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 5)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot6 <- ggplot(df_quantities_part6, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 6)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot7 <- ggplot(df_quantities_part7, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 7)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot8 <- ggplot(df_quantities_part8, aes(x = as.factor(quantities), y = frequency)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 8)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(plot1)
print(plot2)
print(plot3)
print(plot4)
print(plot5)
print(plot6)
print(plot7)
print(plot8)
prescriptions_by_year <- claims_demo %>%
group_by(year) %>%
summarise(count = n())
#1. Bar Graph of Counts by Year
# Plot the bar graph
ggplot(prescriptions_by_year, aes(x = year, y = count)) +
geom_bar(stat = "identity") +
labs(title = "Number of Prescriptions by Year",
x = "Year",
y = "Number of Prescriptions") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
From 2006 to 2013, there was a steady increase in the number of prescriptions, peaking in 2013 and 2014 with over 300,000 prescriptions each year.
This rise suggests a growing trend in prescription rates or improved access to healthcare during these years.
However, following 2014, there is a notable decline in the number of prescriptions, reaching its lowest point around 2019. This decline may be attributed to changes in healthcare policies.
In the subsequent years, there is a slight recovery in 2021 and 2022, indicating a resurgence in prescription rates, but this is followed by a sharp decrease in 2023. This could be influenced by external factors such as changes in healthcare regulations, economic conditions, or the impact of global events like the COVID-19 pandemic.
Filter the data to include only entries with sex “F” and “M”
fig1_filtered <- claims_demo %>%
filter(der_sex %in% c("F", "M"))
# Summarize the data by year and sex
prescriptions_by_year_sex <- fig1_filtered %>%
group_by(year, der_sex) %>%
summarise(count = n())
`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.
# Plot the bar graph
ggplot(prescriptions_by_year_sex, aes(x = year, y = count, fill = der_sex)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Number of Prescriptions by Year, Stratified by Sex",
x = "Year",
y = "Number of Prescriptions",
fill = "Sex") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
The second bar graph shows the number of prescriptions by year, stratified by sex (F and M). It highlights several important trends and differences between male and female prescription rates. Overall, the number of prescriptions for both sexes follows a similar pattern, with an increase from 2006, peaking around 2013 and 2014, followed by a decline.
Throughout the years, males consistently have higher prescription counts compared to females. The peak years of 2013 and 2014 show the most significant differences, with males having substantially more prescriptions than females. After 2014, the number of prescriptions for both sexes declines, reaching a low point around 2019.
This graph indicates that while the overall prescription trends are similar for both males and females, males consistently receive more prescriptions than females across all years.
create a map of prescriptions by US state
fig2<-claims_demo %>% select(pat_state, year)
write.csv(fig2, "fig2.csv")
prescriptions_by_state <- fig2 %>%
group_by(pat_state) %>%
summarise(total_prescriptions = n())
head(prescriptions_by_state)
prescriptions_by_state <- as.data.frame(prescriptions_by_state)
prescriptions_by_state<-prescriptions_by_state[-1,]
colnames(prescriptions_by_state)[1]<-"state"
plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
scale_fill_continuous(name = "Total Prescriptions", label = scales::comma) +
labs(title = "Total Number of Prescriptions by US State") +
theme(legend.position = "right")
state_data <- merge(prescriptions_by_state, state_centers, by.x = "state", by.y = "state")
# Total prescriptions plot using plotly
total_prescriptions_plot <- plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
geom_text(data = state_data, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
scale_fill_gradientn(name = "Total Prescriptions", colors = c("blue", "green", "yellow", "red"), labels = scales::comma,
breaks = seq(0, max(prescriptions_by_state$total_prescriptions), length.out = 6)) +
labs(title = "Total Number of Prescriptions by US State") +
theme(legend.position = "right")
# Convert to plotly object
ggplotly(total_prescriptions_plot)
ui <- fluidPage(
titlePanel("Number of Prescriptions by US State"),
tabsetPanel(
tabPanel("Total",
plotOutput("staticMap")
),
tabPanel("Yearly",
sidebarLayout(
sidebarPanel(
selectInput("year", "Select Year:", choices = unique(prescriptions_by_state_year$year))
),
mainPanel(
plotlyOutput("usMap")
)
)
),
tabPanel("Product Trends",
plotOutput("productTrendGraph")
)
)
)
server <- function(input, output) {
output$staticMap <- renderPlot({
state_data <- merge(prescriptions_by_state, state_centers, by.x = "state", by.y = "state")
plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
geom_text(data = state_data, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
scale_fill_gradientn(name = "Total Prescriptions", colors = c("blue", "green", "yellow", "red"), labels = scales::comma,
breaks = seq(0, max(prescriptions_by_state$total_prescriptions), length.out = 6)) +
labs(title = "Total Number of Prescriptions by US State") +
theme(legend.position = "right")
})
output$usMap <- renderPlotly({
data_year <- prescriptions_by_state_year %>%
filter(year == input$year)
state_data_year <- merge(data_year, state_centers, by = "state")
p <- plot_usmap(data = data_year, values = "total_prescriptions", color = "white") +
geom_text(data = state_data_year, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
scale_fill_gradientn(name = paste("Prescriptions in", input$year), colors = c("blue", "green", "yellow", "red"), labels = scales::comma,
breaks = seq(0, max(data_year$total_prescriptions), length.out = 6)) +
labs(title = paste("Number of Buprenorphine OUD Prescriptions by US State in", input$year)) +
theme(legend.position = "right")
ggplotly(p)
})
output$productTrendGraph <- renderPlot({
ggplot(product_data, aes(x = year, y = percentage, color = product_name)) +
geom_line(size = 1.2) +
scale_color_viridis(discrete = TRUE) +
labs(title = "Percentage of Each Product Over Time (2006-2023)",
x = "Year",
y = "Percentage of Total Prescriptions",
color = "Product") +
theme_minimal() +
theme(legend.position = "right", legend.title = element_text(size = 10),
legend.text = element_text(size = 8))
})
}
shinyApp(ui = ui, server = server)
product_trend_plot <- ggplot(product_data, aes(x = year, y = percentage, color = product_name)) +
geom_line(size = 1.2) +
scale_color_brewer(palette = "Set3") +
labs(title = "Percentage of Each Product Over Time (2006-2023)",
x = "Year",
y = "Percentage of Total Prescriptions",
color = "Product") +
theme_minimal() +
theme(legend.position = "right", legend.title = element_text(size = 10),
legend.text = element_text(size = 8))
# Convert to plotly object
ggplotly(product_trend_plot)
240923*0.572 ## discussion of product Trends 240923-137808 100-57.2 866904/240923
SUBOXONE Trend: There is a significant increase in the percentage of SUBOXONE prescriptions starting around 2010, reaching a peak around 2015-2016, followed by a sharp decline.
Literature Support: SUBOXONE, a combination of buprenorphine and naloxone, was approved by the FDA in 2002 for the treatment of opioid dependence. The increase in its prescription aligns with its effectiveness in treating opioid use disorder (OUD) and its expanded use in clinical practice. The decline in recent years could be attributed to the introduction of new buprenorphine products and the implementation of generic versions, which tend to be less expensive and more accessible.
BUPRENORPHINE HCL Trend: Shows a significant increase starting around 2018. Literature Support: Buprenorphine HCL has been used for both pain management and OUD. The increase in recent years can be attributed to the broader acceptance and understanding of buprenorphine’s efficacy in treating OUD. It may also reflect the FDA’s efforts to expand access to buprenorphine as part of the response to the opioid crisis.
BUPRENORPHINE HCL/NALOXONE Trend: Steady increase in percentage over time after 2013. Literature Support: This trend supports the literature that emphasizes the effectiveness of combination therapies (buprenorphine/naloxone) in reducing the potential for misuse compared to buprenorphine alone. The increase in prescriptions also coincides with the FDA’s encouragement of such combinations for safety reasons.
SUBUTEX Trend: High percentage initially, followed by a decline around 2010 and remaining low thereafter. Literature Support: SUBUTEX, which contains only buprenorphine, was commonly used before the introduction of combination products. Its decline is supported by FDA recommendations favoring combination products like SUBOXONE to reduce misuse potential.
SUBLOCADE Trend: Emerges around 2018 with a rising trend. Literature Support: SUBLOCADE, an extended-release injectable form of buprenorphine, was approved by the FDA in 2017. Its increasing prescription rate reflects the growing acceptance and utilization of long-acting formulations for better compliance and sustained recovery in OUD treatment.
PROBUPHINE IMPLANT KIT Trend: Shows a very low and relatively flat trend. Literature Support: PROBUPHINE, approved in 2016, is an implantable form of buprenorphine designed to provide continuous delivery over six months. Its limited use might be due to the invasive nature of implantation and the preference for less invasive administration routes.
ZUBSOLV Trend: Shows a minor presence with some increase over time. Literature Support: ZUBSOLV is another buprenorphine/naloxone combination product with a different formulation and bioavailability. Its introduction provided an alternative to SUBOXONE, and I believe is highly reliant on prescriber habits.
BUNAVAIL Trend: Minor presence with slight fluctuation. Literature Support: BUNAVAIL, a buccal film formulation of buprenorphine/naloxone, was introduced as an alternative to sublingual formulations. Its relatively minor uptake could be due to competition with other established combination products.
Trends Not Fully Supported by Literature
Sharp Decline in SUBOXONE: The rapid decline in SUBOXONE prescriptions post-2015 is sharper than might be expected solely due to the introduction of generics and new formulations. This could indicate other factors at play, such as changes in insurance coverage, prescriber preferences, or regulatory actions targeting specific products.
Relatively Flat Trend for PROBUPHINE: Despite its innovative implant approach, PROBUPHINE’s limited adoption might not fully align with its potential benefits as highlighted in clinical studies. Factors such as cost, patient and provider acceptance, and procedural complexity likely play a role.
rectype_summary <- summary(claims_demo$rectype)
dosage_form_nm_summary <- summary(claims_demo$dosage_form_nm)
route_summary <- summary(claims_demo$route)
strength_summary <- summary(claims_demo$strength)
rectype_plot <- ggplot(claims_demo, aes(x = factor(rectype))) +
geom_bar(fill = 'skyblue', color = 'black') +
labs(title = 'Distribution of Rectype', x = 'Rectype', y = 'Frequency') +
theme_minimal()
dosage_form_nm_plot <- ggplot(claims_demo, aes(x = factor(dosage_form_nm))) +
geom_bar(fill = 'salmon', color = 'black') +
labs(title = 'Distribution of Dosage Form Name', x = 'Dosage Form Name', y = 'Frequency') +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
route_plot <- ggplot(claims_demo, aes(x = factor(route))) +
geom_bar(fill = 'lightgreen', color = 'black') +
labs(title = 'Distribution of Route', x = 'Route', y = 'Frequency') +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
strength_plot <- ggplot(claims_demo, aes(x = factor(strength))) +
geom_bar(fill = 'lightblue', color = 'black') +
labs(title = 'Distribution of Strength', x = 'Strength', y = 'Frequency') +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
rectype_plot
dosage_form_nm_plot
route_plot
strength_plot
colnames(claims_demo)
selected_data <- claims_demo %>% select(pat_id, der_sex, pat_region, pat_state, yob_group, zip3_group)
Error in `select()`:
! Can't subset columns that don't exist.
✖ Column `yob_group` doesn't exist.
Backtrace:
1. claims_demo %>% ...
3. dplyr:::select.data.frame(., pat_id, der_sex, pat_region, pat_state, yob_group, zip3_group)